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We present a dynamical model of DNA mechanical unzipping under the action of a force. The 
model includes the motion of the fork in the sequence-dependent landscape, the trap(s) acting on 
the bead(s), and the polymeric components of the molecular construction (unzipped single strands 
of DNA, and linkers). Different setups are considered to test the model, and the outcome of the 
simulations is compared to simpler dynamical models existing in the literature where polymers are 
assumed to be at equilibrium. 




I. INTRODUCTION 



Over the past fifteen years various single molecule ex- 
periments have investigated DNA mechanical and struc- 
tural properties Jl|,J| S i, S, S 0, i, i, 0, HH, E 

rotein-DNA interactions 
23, m, m. These exper- 
iments provide dynamical information usually hidden in 
large scale bulk experiments, such as fluctuations on the 
scale of the individual molecule. The separation of the 
two strands of a DNA molecule under a mechanical stress, 
usually referred to as unzipping, has first been carried out 
by Bockelmann and Heslot in 1997 [Sj. The strands are 
pulled apart at a constant velocity while the force nec- 
essary to the opening is measured. The average opening 
force for the A phage sequence is of about 15 pN (at 
room temperature and standard ionic conditions), with 
fluctuations around this value that depend on the partic- 
ular sequence content. Bockelmann, Heslot and collabo- 
rators have shown that the force signal is correlated to 
the average sequence on the scale of ten base-pairs but 
could be affected by the mutation of one base-pair ad- 
equately located along the sequence figl . Liphart et al. 
and [l^ and Danilowicz et al. have performed 

an analogous experiment, using a constant force setup, 
on a short RNA and a long DNA molecules respectively 
(Figure [TJ3). The distance between the two strands ex- 
tremities is measured as a function of the time while the 
molecule is submitted to a constant force. DNA strands 
separation has been also studied in single molecule ex- 
periments by translocation through nanopores (26l . Wf\ . 

The motivation underlying unzipping experiments of 
DNA is (at least) two-fold. First the study of unzip- 
ping aims at a better understanding of the mechanisms 
governing the opening of DNA during transcription and 
replication by proteins such as Polymerases, Helicases, 
and Exonucleases [l^, [2ll, [2^ [2^. Simple theoretical 
models describing the opening as an unidimensional ran- 
dom walk on a sequence-dependent free energy landscape 
have been proved to describe quite well several experi- 
mental effect such as stick-slip motion in the opening at 



constant velocity [T3| , the lon g pa uses at fixed position 
of unzipping at constant forces [l6l [sol HH , the hopping 
dynamics between two or more states in unzippin g a t 
critical forces of short DNA molecules [H, [Sll, [H, 111] , 
and torsional drag effects in unzipping at large velocity 
[Til [s^l- Moreover statistical mechanical analysis have 
been successfully applied to extract from experimental 
data the sequence-dependent free energy landscape and 
the height of free energy barriers [H, [36| . 

Secondly unzipping experiments could potentially be use- 
ful to extract information on the sequence itself [sj- 
Recently single molecule sequencing has been achieved 
by monitoring a DNA/RNA polymerase in the course of 
DNA synthesis from a ssDNA template [H, ; such sin- 
gle molecule sequencing could become competitive with 
standard DNA sequencing because they do not require, a 
priori, amplification through polymerase chain reactions. 
A fundamental question on the possibility of extracting 
information on the sequence from unzipping experiments 
is the influence of the experimental setup on the mea- 
sures and the limitations imposed by the latter [s^, HI] . 
Indeed characteristic spatio-temporal limitations are the 
flnite rates of data acquisition, the relaxation time of the 
bead, the limited spatial resolution, the thermal drift, 
and more generally the noise in the instruments. More- 
over the dynamics of the opening fork (Figure [T]) is influ- 
enced by the single strands (open parts) of the molecule 
and the linkers, and cannot be deduced directly from the 
observation of the bead from which the force or the po- 
sition are measured. 

The accuracy of unzipping experiments at flxed velocity 
have improved a lot over the last decade. Initially per- 
formed with an optical fiber Q, experiments were then 
based on the use of sim ple optical traps [13] ■ Nowadays 
double optical traps [l3l.l36| allow to reduce considerably 
the drift of the setup and to achieve a temporal reso- 
lution of the order of 10 KHz, a sub-nanometric spatial 
resolution, and a precision on measured forces of the or- 
der of fraction of pN. Unzipping at fixed force has been 
performed by a magnetic trap with a low temporal reso- 
lution (from 60 Hz to 200 Hz) due to the time needed to 



extract the position of the bead, the spatial precision be- 
ing of the order of 10 nm/-\/Hz) (28l.[29l|. or by an optical 
trap also with a low temporal resolution (about 10 Hz) 
imposed by a feedback mechanism needed to keep the 
force constant |15|. Recently a new dumbbell dual opti- 
cal trap has been developed. It operates without feed- 
back and can maintain the force constant over distances 
of about 50 nm [s^ with temporal resolution of 10 kHz 
and a spatial resolution of 0.1 nm/-\/Hz. 

Limitations due to the experimental systems were first 
addressed in fsoj. This paper stated the impossibility to 
infer the sequence due to ssDNA fluctuations: fluctua- 
tions increase with the number of opened base-pairs and 
can become larger than the length of about 1 nm corre- 
sponding to the spatial resolution of one open base-pair. 
This problem could however be solved by integrating out 
the single strand dynamical fluctuations. Several works 
have studied the effects of the setup on the hopping dy- 
namics of small RNA molecules [sl, IH, H^, |40]. The 
following effects have been underlined. First of all the 
free energy landscape changes when adding to the free 
energ y an harmonic potential, due to the bead and han- 
dles, [Tol [3^ [3^ H^l- Therefore, for a given force, the 
measured separation of the extremities depends on the 
stiffnesses of the trap and handles. Moreover, the open- 
ing and closing rates depend on the stiffness of the opti- 
cal trap; in particular when the experimental system gets 
softer the fluctuations of the force gets smaller, and the 
hopping rates approach their fixed- force values. 

In this paper we introduce a model for the coupled dy- 
namics of the opening fork, the ssDNA strand, the link- 
ers, and the bead in the optical or magnetic trap. Essen- 
tial notions and existing literature are reviewed in Section 
II. Our dynamical model is presented in Section III. Our 
program allows to simulate a generic setup, character- 
ized by bead dimensions, optical stiffness (absent in case 
of magnetical tweezers), linker composition (dsDNA or 
ssDNA) and lengths, length of molecule to be unzipped. 
All the parameters that characterize the different dynam- 
ical components can be adjusted in the simulation. The 
model is then used to simulate fixed-force (Section IV) 
and fixed-extension (Section V) numerical unzippings. 
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FIG. 1: Typical experimental setups that will be described 
in the following. A) A setup with two optical traps (beads 
xi and X4) drawn as springs and whose centers are the black 
vertical lines; B) a setup with a single magnetic bead X3 that 
applies a constant force on the molecule attached to a fixed 
"wall". In both cases the molecular construction is made by 
a DNA molecule that has to be opened (therefore one should 
include two single-strand linkers that are the opened parts 
of the molecule) and one double-stranded DNA linker. The 
coordinates Xi are the distances of the corresponding points 
from the left reference position (which is the center of the left 
optical trap in case A and the fixed wall in case B). 



A. Thermodynamics of the components 

1. Polymeric models for the linkers and open molecule 

A polymer model is specified by its free energy as a func- 
tion of the extension x for a given number n of monomers; 
we call this quantity W{x,n). When x and n are large 
W is an extensive quantity hence W{x,n) = nw{x/n) = 
nw{l), where I = x/n is the extension per monomer. We 
also define 



lif) 



dW{x,n) 



w'{l) 



dx 

inverse of f{l) , 



(1) 



g{f) = max[/ 1 - w{l)] ^ f l{f) - w[l{f)] 



II. FREE-ENERGIES, TIME SCALES, AND 
EFFECTIVE DYNAMICS 



We discuss hereafter the thermodynamic properties of 
the various parts of the experimental set-up (DNA se- 
quence, open part of the molecule, single- or double- 
strand linkers), as well as the relevant time scales. Finally 
we briefly review previous dynamical studies where the 
linkers and the open portion of the molecules are assumed 
to be at equilibrium. 



which are, respectively, the force at fixed extension, the 
average extension at fixed force, and the free-energy at 
fixed force. Notice that g{f) is simply the integral of 
1(f)- Hence a polymer model is completely described 
from the knowledge of the extension versus force char- 
acteristic curve, l{f). In the following we will use some 
classical models for this function: 



Gaussian (Hook) model: 



I'Hook 
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TABLE I: Binding free energies <]o(bi, bj+i) (units of fcsT) 
obtained from the MFOLD server [13, El] for DNA at room 
temperature, pH=7.5, and ionic concentration of 0.15 M. The 
base values bi,bi+i are given by the line and column respec- 
tively. 



where the stiffness constant fc™ is related to the 
temperature T and the average squared monomer 
length (at zero force) through fc™ — ksT/b^. 

• Freely-Jointed Chain (FJC) model: 

fb \ ksT 



Ipjcif) = coth 



fb 



(3) 



is the extension (per monomer) of a chain of rigid 
rods of length 6, free to rotate around each other. 
Comparison of this model with force-extension 
curves for single-stranded DNA shows that a better 
fit is obtained from a Modified FJC, 



lMFJc{f)^d(l + 



f 



Is 



X I 



FJC 



if) 



(4) 



which takes into account elasticity effects on the 
rod length. Standard fit parameters are d = 0.56 
nm, b = 1.4 nm, and jss — 800 pN. 



Extensible Worm-Like Chain (WLC) model: 



IwLcif) — L 



1 

1 - - 
2 



keT 
fA 



1/2 



Ids 



(5) 



is the formula for the high-force extension of an 
elastic chain with persistence length equal to A. 
Experiments shows it is an excellent description of 
double-stranded DNA at high forces, with L = 0.34 
nm, A = 48 nm and -yds = 1000 pN. 



2. Free-energy landscape for the sequence 



Let bi = A, T, C, or G denote the i**^ base along the 
5' — > 3' strand (the other strand is complementary), and 
B = {bi,b2, ■ ■ ■ ,b]y}. The free energy excess when the 
first n bp of the molecule are open with respect to the 
closed configuration (n = 0) is [31| 



G{n;B) ^^ga{bi,bi+i) 



(6) 
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FIG. 2: Free energy G (units of k_BT) to open the first n base- 
pairs, for the first 50 bases of the DNA A-phage at forces 15.9 
(dashed curve) and 16.4 pN (full curve). For / = 15.9 pN the 
two minima at bp 1 and bp 50 are separated by a barrier of 
12 kgT. Inset: additional barrier representing the dynamical 
rates (|2ip to go from base 10 to 9 (barrier equal to 2gss='2-5 
ksT), and from base 9 to 10 (barrier equal to (;o(&9, 6io)=3 
kflT), see text. 



go{bi,bi+i) is the binding energy of base-pair (bp) num- 
ber «; it depends on bi (pairing interactions) and on the 
neighboring bp due to stackirig interactions. is 
obtained from the MFOLD server [4l|, |42| , and listed in 
Table U for a 150 mM NaCl, room temperature and pH 
7.5. The values of the free-energies should be changed 
for different ionic conditions and temperatures. 



As an illustration we plot the free energy G{n\ A) of 
the first 50 bases of the A phage sequence, A = 
(Ai, A2, . . . , Ajv), in Figure Rafter subtraction of n gss{f) 
for forces / — 15.9 and 16.4 pN. gss{f) is the work to 
stretch the two opened single strands when one more bp 
is opened, and calculated from the Modified FJC model 
(|3]). The substraction allows to compare the increase in 
the free-energy due to the opening of the sequence to the 
gain resulting from the release of ssDNA polymers at a 
given force. 



At these forces the two global minima in Figure [2] are lo- 
cated in n = 1 (closed state) and n = 50 (partially open 
state). Experiments on a small RNA molecule, called 
P5ab have been performed at the critical force fc 
such that the closed state has the same free energy than 
the open one: G{N;A) = N gssifc)- They showed that, 
as the barrier between these two minima is not too high, 
the molecule switches between these two states, see Sec- 
tion HTCl 
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TABLE II: Fluctuations of single stranded DNA at / = 15 
pN and T = 16.7 °C; 5x/x = 0.37/ ^n, Sf/f = 1.57/Vn, 
T = 4.83 • 10"" s 
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TABLE III: Fluctuations of double stranded DNA at / = 15 
pN and T = 16.7 °C; 5x/x = 0.17/ Vn, 5f/f = 4.83/Vn, 
T = 5.1 ■ 10"^^ s 



B. Fluctuations at equilibrium 

1. Case of a single polymer 

We now consider the orders of magnitude of the fluc- 
tuations of the polymer. When submitted to a force 
of / = 15 pN, the average extension of the polymer is 
X = nx™ with X™ = l{f). We use for single-stranded 
DNA the MFJC model, and for double-stranded DNA 
the WLC model, with the parameters discussed in sec- 
tion lll ATI then we get = 0.46 nm and xJJ^ = 0.33 nm 
for the ss- and ds-DNA. At thermal equilibrium the ex- 
tension will fluctuate around these average values. The 
fluctuations are controlled by the microscopic effective 
spring constant k''"'{l) = w"{l) — l/l'{f). For ds- and 
ss-DNA we find, respectively, fc^ — 1311 pN/nm and 
fc™ =138 pN/nm according to the models above. For 
a polymer with n monomers the stiffness is fc = k"^/n 
since the effective spring constant is given by k{x,n) — 
£,W{x,n) = k'^{x/n)ln. 

Alternatively the force / exterted on the polymer will 
fluctuate around its average value / if its extremities are 
kept at a fixed distance x from each other. These fluc- 
tuations of force (in the fixed-extension ensemble) and 
extension (in the fixed-force ensemble) are easily com- 
puted by a quadratic expansion of the free-energy around 
the average i.e. when approximating the polymer with a 
spring of stiffness k™/n, with the result 



(Sx'^) = —5 — n 



kgT k''' 



(7) 



Defining Sx — y/J^jp) and 5f = \/(5p) we get 



5x 



knT 



f 



knTk"" 1 



P 



(8) 



As expected the relative fluctuations of both force and 
extension become smaller and smaller as the number n of 
monomers increases. Some values are reported in Tables 

iniandiini 



2. Case of several polymers (fixed- distance setup) 

Consider now the case of several polymers e.g. linker and 
open part of the molecule attached one after the other. In 
a fixed-force experiment the components of the setup are 
independent (at the level of the saddle-point approxima- 
tion), and the fluctuations in the extensions simply add 
up. In the fixed-distance setup, however, correlations be- 
tween the extension make the analysis more complicated. 
As a concrete example we consider the setup in fig- 
ure[TJ\. The linker joining xi and X2 is a double-stranded 
DNA segment of Nds bases. The two linkers joining 
{x2,X3) and {x3,X4) are single-stranded DNA segments 
of Nss = N^s "I" bases, where n is the number of 
opened base-pairs. The centers of the two optical traps 
are in and X. We call xi the position of the first 
bead, and X4 the position of the second. The probability 
Peq{n, xi, X2, X3, X4) = e~^/'''^'^ where the free energy F 
reads 

F{x,n) = ^kixl + Wds{x2 - xi,Nds) 

+ Wss{x3 -X2,Nss) + Wss{x4 - X3, Nss) (9) 

+ h2ixi-Xf + G{72;B) . 

where Wds{x,Nds) = NdsWds{x/Nds) and Wss{x,Nss) = 
NssWss{x/Nss) are the elongation free energies of the 
double strand and single strand respectively. 
In order to study the fluctuations in this setup we first 
find the maximum of Peq assuming that G(n; B) — n qq 
i.e. a uniform sequence B, and treating n as a continuous 
variable assuming that it is large. At the maximum Xi — 
Xi and we define 



X2 - Xi 

Nds 

X3 - X2 
Nss 



Xj - X3 
Nss 



(10) 



The saddle-point condition dxiFA = gives the follow- 
ing equations, that represent the force balance condition 
along the chain: 

fclil = W'dsix^) = w'ssixTs) ^k2iX- X,) EE f . (11) 

The derivative with respect to n gives, using Eqs. (jlip 
and ll]), the condition 



go = 2[xTsw'ss{x:i) ~ WssixTs)] = gUf) 



(12) 
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TABLE IV: Saddle-point calculation for the setup in figure[TJ\. 
with a uniform molecule and fci = 0.1 pN/nm, k2 = 0.512 
pN/nm, Nds = 3120, iVf^ = 40. The force along the molecule 
is / = 16.5, then fc™ = 152 pN/nm, fcj^ = 1416 pN/nm and 
fclft = 0.07 pN/nm. 



that allows to find the force / transmitted along the 
chain. Once is solved, the extension of the beads 
and of the double and single stranded part of DNA (xi, 
AT — 54, and respectively) are determined by 
Eq. (fTT|) . Finally the number of open bases n is deter- 
mined by 

xi + Ndsx"al + 2(7V°, + n)xZ + {X - xa) ^ X . (13) 

Note that the value of / is determined only by go- 
We work at temperature T = 16.7 °C (fc^T = 4 pN 
nm) and choose a uniform molecule with go — 2.69 fesT, 
which is a representative value for the pairing free ener- 
gies in tableU We use the same models as in section lTlB II 
for the single and double stranded DNA, with Nds = 3120 
and iV°, = 40. Then solving Eq.dn]) we get / = 16.5 pN, 
and from Eq. pTj) we get = 0.47 nm, x^^ = 0.33 nm. 
We choose fci = 0.1 pN/nm, then xi = 165 nm, and 
k2 = 0.512 pN/nm, then AT — 5:4 = 32 nm. Given these 
values, n is defined by X using Eq. 



X - 1264 



(14) 



with X expressed in nanometers. 

For the same setup we can compute the fiuctuations of 
n and of the elongations of the elements of the setup. 
In particular the fluctuations of the bead positions are 
measurable in the experiment. 

Let us define Sxi — Xi — Xi and Sn ^ n ~ n. To simplify 
the formalism we also define Sxds = 5x2 — 6x1, Sx^^ = 
6x3 — 5x2, and Sxf^ = 6x4 — 6x3. A quadratic expansion 
of F around its minimum gives 



SF^ 



^kM + lk25xl^^^^^^ 



5x1 



2N° 4 



-[i5xi 



2Nas 
xTJnf + i5xf. 



(15) 



I6n) 



Using (HI) and ^ we get = w'^,{x2) = 152 pN/nm 
and kZ = w'^,{xZ) = 1416 pN/nm. 
One should take care of the fact that 5xi + Sx^ + 5xds + 
dXgg + 5xfg = 0; it is convenient to express dxf^ as a 
function of the others since its fluctuations are identi- 
cal to the ones of Sx^^. The quadratic expansion of 
the function 5F has the form 5F = ^SxASx where 



(5x = {5xi, 5x4, 5xds, 5x 
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The inverse of the matrix A is 
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This immediately gives 
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(16) 



(17) 



(18) 



(19) 



4fceff(x™)2 

and shows that the fluctuations of n are dominated by 
the weakest element of the setup; moreover the corre- 
lation between the beads displacements 5xi, 5x4 and 
the fluctuations of the number of open base-pairs 5n is 
{5n5xi) = -2!^ and (^^4) = -^^j the stiffer 
the optical trap , the weaker is the correlation between 
the location of the bead and the number of open bases. 
Examples are given in Table ITVl 



C. Effective dynamical models 

In the simplest dynamical models the fork (separating 
the open and closed portions of the molecule) undergoes 



6 



a biased random motion in the sequence landscape. The 
hnkers are treated at equilibrium, which is correct if their 
characteristic time scales are much smaller than the av- 
erage time to open or close a base-pair. 



1. Time scales for the polymeric components of the setup 

In this section we recall the typical time scales of the 
polymeric components in the setup. Assume that the 
polymers are subject to a Brownian force ri{t) which is a 
zero-average Gaussian process with autocorrelation func- 
tion {r]{t)T]{0)) = 2TT5{t). Let T be the friction coeffi- 
cient of the polymer ^43|], that is, the ratio between the 
viscous force exerted by the solvent and the velocity. As 
will be shown in Section UTTl the friction coefficient scales 
as r = 7™n/3 with 7™ = 7^^ - 2-10-8 pN s/nm^. Then, 
approximating f{x,n) ~ k''^x/n, the relaxation time for 
an isolated polymer of n bases is given by 
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FIG. 3: Number of open base-pairs as a function of the time 
for various forces (shown on Figure). Data show one numer- 
ical unzipping (for each force) obtained from a Monte Carlo 
simulation of the random walk motion of the fork with rates 



7"'n 



(20) 



Note that the factor 3 at the denominator of the above 
equation is an approximation for the true factor 7r^/4. 
The validity of its approximation, and the simplification 
it leads to will be discussed in Appendix VK\ 
It is useful to compare the amplitude of the force fluc- 
tuations with the noise. To do this we approximate 
{Sf{t)Sf{0)) ^ 2T{Sp)Sit) = 2TTfS{t). Then, using 
Eq. dH) to estimate {6f), we get F/ = n7™/3 = F, and 
(not surprisingly) the force fluctuations are of the same 
order of the noise term. 

From TablelUthe relaxation time of the unzipped strands 
is smaller that the typical base-pair opening (or closing) 
time as long as the number n of unzipped bases is smaller 
than a few hundreds. This is the case in particular for 
unzipping experiments on short RNA molecules. 



2. Random walk in the sequence landscape 

Let us first model the motion of the fork alone, that 
is, assuming that the other components of the setup 
are at equilibrium. We consider a DNA molecule un- 
zipped under a fixed force / in the sequence-landscape 
G{n; B) — n gssif) of Figure O The fork, whose position 
is denoted by n{t), can move forward (n —^ n + 1) or 
backward (n — > n — 1) with rates (probability per unit of 
time) equal to, respectively. 



ro(6„+i, 6,1+2) = r exp [ - l3go{bn+i,bn+2)] , 
rc = r exp [-2Pgss{f)] , 



(21) 



where f3 = l/fc^T, see Fig[21 The value of the attempt 
frequency r is of the order of 10^ Hz fH, [H [3l|. Ex- 
pression (PT|) for the rates is derived from the following 



assumptions. First the rates should satisfy detailed bal- 
ance. Secondly we impose that the opening rate Tq de- 
pends on the binding free energy, and not on the force, 
and vice- versa for the closing rate rc- This choice is mo- 
tivated by the fact that the range for base-pair interac- 
tion is very small: the hydrogen and stacking bonds are 
broken when the bases are kept apart at a fraction of 
Angstrom, while the force work is appreciable on the dis- 
tance of the opened bases (« 1 nm). On the contrary, 
to close the base pairs, one has first to work against the 
applied force, therefore the closing rate rc depends on the 
force but not on the sequence. This physical origin of the 
rates is reported in the inset of Figure [51 Notice that, as 
room temperature is much smaller than the thermal de- 
naturation temperature, we safely discard the existence 
of denatured bubble in the zipped DNA portion. 

An example of unzipping dynamics for the A-phage se- 
quence is shown in Figure [3] The characteristic pauses in 
the unzipping, present in experiments and corresponding 
to deep local minima in the sequence-landscape are re- 
produced. The rates (j21|) lead to a master equation for 
the probability Pn(t) for the fork to be at site n at time 
t 



dt 



N 



(22) 



where the matrix r„.m is tridiagonal with nonzero en- 
tries Trn-i.m = -rdf), Tm+i,™ = ~ro{in) and T^.m = 
7'o(m,) -f- rdf). Given this transition matrix the open- 
ing dynamics can be simulated with a Monte Garlo. For 
small RNA or DNA molecules the transition matrix T„.,„ 
can be diagonalized numerically [sij. The smallest non- 
zero eigenvalue gives the switching time between a closed 
and open configuration for hairpin with a free energy bar- 
rier such as the one plotted in Figure [H 
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Dynamics of the bead with equilibrated linkers ana 
strands 



In a typical experiment the force is exerted on the 
molecule through the action of a (magnetic or optical) 
trap on the bead. While the external force on the 
bead can be considered as constant (e.g. in a mag- 
netic trap), the force acting on the fork fluctuates un- 
less the trap (and the molecular construction) is very 
soft, see Eq. ([5]). Therefore the fixed- force model of the 
previous section has to be modified. In addition the 
bead, of size i? ~ 1 ^m, is a slow component whose 
dynamics need to be taken into account. Let us de- 
note by k the stiffness of the trap and 7 the friction 
of the bead in the solvent of viscosity 77. Typical val- 
ues for these quantities are k — 0.1 — 0.5 pN/nm, and 
7 = QirRrj = 1.67 10^^ pN s/nm. The characteristic 
relaxation time of the bead is thus t — ^/k ^ 0.2 — 1 ms. 
The coupled dynamics of the fork and the bead was con- 
sidered by Manosas et al. [ij] in the case of small RNA 
unzipping, with a single optical trap. For such small 
molecules the relaxation time of the unzipped strands is 
expected to be much smaller than the characteristic time 
of the bead, and the molecule can be considered at equi- 
librium. The dynamical scheme therefore consists in a 
coupled evolution equation for the location of the bead 
and of the fork. The bead position obeys a Langevin 
equation including the external force and the force ex- 
erted by the fork through the (equilibrated) linkers and 
unzipped strands, while the fork moves with rates (j21[) 
with a bead location-dependent force. 
A main conclusion of Ref. [l3| is that, in the absence 
of feedback imposing a fixed force on the molecule, the 
trap stiffness must be as low as possible to detect jumps 
between closed and open configurations of the RNA 
molecule. We will discuss the validity of this statement 
in an information-theoretic setting in section fVBI 



III. DYNAMICAL MODELLING OF THE 
SETUP AND ITS COMPONENTS 

The assumption that the linkers and the unzipped 
strands are at equilibrium as the unzipping proceeds is 



correct for short molecules as was the case in [IJ]. For 
long DNA molecules the relaxation time of the unzipped 
strands may become large, and a dynamical modelling 
of the polymers involved in the molecular construction 
cannot be avoided. 

The purpose of this section is to describe how such a dy- 
namical model can be implemented. We hereafter denote 
by "setup" the full molecular construction that is used in 
a given experiment, including linkers, beads, etc. while 
the word "molecule" refers to the part of DNA which 
has to be opened. In an idealized description, the state 
variable is a vector x = (xi, • • • , Xp) whose elements are 
the distances from a reference position (that can be ei- 
ther the center of an optical trap, or a fixed 'wall' to 



which the polymers are attached) of the extremities of 
the polymeric components in the setup. In addition to a?, 
the number of open base-pairs n is needed to complete 
the description of the state of the setup. 
As discussed in section Hi Al thc total free energy F{x,n) 
of a setup is the sum of different contributions coming 
from all the elements of the setup. A typical example is 
given in Eq. 

Our aim is thus to construct a dynamical model that 
holds on intermediate time scales, t > 10"^ s, and 

1. gives the correct equilibrium Gibbs measure 
Peq{x,n) = exp(-i^(x,n)/(fcBT)), 

2. reproduces the relaxation times for the different el- 
ements of the setup, as discussed below, 

3. gives reasonable dynamical correlations between 
different elements of the setup. 

It is worth to stress at this point that ours is a coarse- 
grained model which does not take into account the mo- 
tion of the individual monomers. It is expected that the 
dynamics on time-scales smaller than the typical sojourn 
time of the fork on a base (> 10~® s) is not relevant to 
our study of unzipping. 



A. Langevin dynamics for the polymers and the 
beads 

First we consider the dynamics of x at fixed n. In Ap- 
pendixl^we show that for long enough times the dynam- 
ics of the setup can be described by a system of coupled 
Langevin equations: 



Tij Xj 



OF 
dxi 



(23) 



where i, j = 1, • • • ,p and: 

• the free energy F{x) is the sum of a contribution 
coming from each element of the setup: 

— each optical trap contributes ^kAx^, where 
Aa; is its elongation; 

— a bead in position i subjected to a constant 
force gives a contribution —fxi] 

— a polymer gives a contribution Wi{Ax, Ni), 
with Ax its elongation and iV^ its number of 
monomers. 

For example, the total free energy of the setups in 
figure [1] are 

Fa{x) = ^kixj + Wds{x2 - Xi,Nds) + Wss{x3 - X2,Nss) 

+ Wss{x4 - X3, Nss) + ^k2{x4 ~ Xf , 

Fb{x) = Wds{xi,Nds) + Wss{x2 - Xi,Nss) 
Wss{xs - X2,Nss) - fx3 . 

(24) 
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• 77 is a Gaussian white noise with zero average and 
variance {'r]i{t)Vj{^)) = 2kBTrijS{t), as requested 
by the fluctuation-dissipation relation; 

• the matrix F is a tridiagonal matrix such that: 

— the diagonal element Ta is the sum of three 
contributions: 

* a term 7™ ;^iVi_i/3-|-7™7Vi/3 coming from 
the adjacent polymers (if any); 

* a term 7 coming from the bead (if any) 
attached to x^; 

* a term taking into account the viscosity 
of the Nc base-pairs of the DNA molecule 
attached to the fork (X3 and X2 in fig- 
ure [TJ'V and B respectively) that are not 

open. This term has the Fleury form 

3/5 

7mo; = I'Nc and has to be added to 
the diagonal element of T corresponding 
to the fork position; 

— the off diagonal elements are zero, except 



I'i+i 'q" that get a contri- 



bution from the polymer joining Xi and x 



i+i- 



For instance the setups in figure [T] correspond to 
the matrices: 



' Ids 3 



IBS 6 



Ids 6 



V 



3 



27, 



7; 

m Was 

3 

7r 



3/5 



< 6 





iss 6 

1 + iTs^j 



7d' 



o\ 



(25) 



A detailed derivation of these results and in particular of 
the form of the matrix F can be found in Appendix |^ 



B. Fork dynamics 

The Langevin equation for the polymer dynamics at fixed 
n must be complemented with transition rates for the 
dynamics of n. To this aim we discretize the Langevin 
equation with time step At, and at each time step we 
allow the opening n — > n -I- 1 or closing — > n — 1 of a 
base-pair at most. 

The dynamics takes the form of a discrete time Markov 
chain, with transitions {x,n) —^ {x',n') and n' S {n,n± 
1}. The total free energy F{x,n) = Fsotup(^, 'T') + 
G{n; B), where the first contribution has been discussed 
in the previous section and G{n] B) is the pairing free 
energy of the molecule, as discussed in section III A 21 
In Appendix [B] we show that in order to satisfy the 
detailed balance condition with respect to Peg (x, n) = 



exp{—F{x,n)/{kBT)), one should perform a single step 
following the procedure: 

1. Choose whether to stay (n' = n), to open (n' = 
n + 1) or to close (n' = n — 1) a base, with rates 
r^'°''^{x,n) respectively: 

r°(f,n)=rAte'^['^("^^)-«("+i;^)l , 
r%x, n) = rMe^^^'^'^^-P^^^'^-^^ , (26) 
(a', n) = 1 — r° (x, n) — r'^ (x, n) . 

2. If the choice was to open, first perform a discrete 
Langevin step x — > x' at fixed n and then increase 
n by one. 

3. If the choice was to close, first decrease n by one 
and then perform a discrete Langevin step x — > x' 
at fixed n' = n — 1. 

4. If the choice was to stay, just perform a discrete 
Langevin step x — s- x' at fixed n. 

The Langevin equation is discrctizcd in a standard way 
by integrating Eq. over a time At: 



x,{t + M) =x,(At)+F,r.i 



, (27) 



where Ej = r]j{t)dt are Gaussian variables with zero 
average and variance 

{E,Ej) = 2kBTT,jAt (28) 

that are independently drawn at each discrete time step. 

C. Free energy at finite n 



In section III Al we discussed some models for the free 
energy W{x,n) of a polymer with n monomers and ex- 
tension X. In the limit x,n ^ 00 at fixed extension per 
monomer, I — x/n, the free-energy enjoys an extensivity 
property: W{x, n) = n w(l). However, in our simulations 
we might be interested to regimes where n is small, typi- 
cally of the order of 10 — 40 for small RNA molecules. In 
this case the knowledge of the free-energy per monomer, 
w, is not sufficient, and a more detailed expression for W 
is necessary to avoid inconsistencies. 
As a starting point of the analysis, we consider a polymer 
made of N indentical monomers whose endpoints are de- 
noted by Ui , i = 1 • • • TV with uq = 0. The Hamiltonian of 
the chain is the sum of pairwise interactions ip{ui — Wi_i) 
and the free energy reads, for x — u^: 



-l3W{x,n) _ p-N+\ 

e — Iq 



dui ■ ■ ■ duN-ie 



(29) 
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where £o is a reference microscopic length scale. From 
the above relation it follows the Chapman-Kolmogorov 
equation 



- f3W (x ,n+m) 



*0 



dye 



-f^W{y,n)—l3W(x—y,7n) 



(30) 



We consider first for simplicity the Gaussian model, 



ip{x) = ifc" 

W{x, n) 



Then it is easy to show that 



2n 



log 



2'KkBT n 



(31) 



In the limit of large polymers one obtains the free-energy 
of a monomer of extension I through 



1 



'w{l) = lim — W{x — ln,n) — (p{l) 



(32) 



as expected and consistently with the discussion of sec- 
tion ITTX] The logarithmic term in (|3ip does contribute 
neither to w, nor to the Langevin equation for x. How- 
ever it does contribute to the rate to close a base-pair, 
see Eq. and should be taken into account in order 

to recover the correct rates. An example of the effect 
of this term is obtained by computing the equilibrium 
probability of n. Consider the (unrealistic) case of a ho- 
mopolymer, G{n; B) = ngo, subject to a constant force 
and using a Gaussian model for the open part of the 
molecule; then 



P. 



(33) 



Therefore Peqiji) is a pure exponential, while if the cor- 
rection were neglected one would have obtained a wrong 
behavior at small n. 

For a generic model of ip{x), one cannot compute 
W{x,n). Still we found that for our purposes (n > 40) 
a consistent approximation is obtained by keeping only 
the first correction to the n — > cxd result, i.e. by defining: 



-l3W{x,n) _ -0nw(x/n) Pk{x / n)l1 



27112 



(34) 



where k{l) = w"{l). One can check that this expres- 
sion satisfies Eq. ([50]) with corrections in the exponent 
of 0(1), while the terms 0(logn -I- logm) are taken 
into account. Within this approximation, the error in 
logrc(x,n) in Eq. (pS)) is 0(l/n^) while if the first cor- 
rections are neglected it is 0{l/n). 

In the following we will make use of the definition (p4| un- 
less otherwise stated. We will discuss an example where 
the effects of neglecting the corrections is clearly observ- 
able. 



The total free energies of the two setups are given 
by Eq. ^ plus the term G(n; B). 



• The free energy of each polymer includes the 
saddle-point corrections, i.e. it is given by Eq. (IM)) . 
The relation l{f), see section III A| is numerically 
inverted to obtain wll) and k{l) that enter in 
Eq. (121 • 



• For the single-stranded DNA we used the MFJC 
model, Eq. ([4]), with d = 0.56 nm, b — lA nm and 
7,, = 800 pN. 

• For the double-stranded DNA we used the WLC 
model in Eq. ([5|), with a small regularization term 
to avoid a divergence for / — )■ 0, which is however 
irrelevant for values of forces to be discussed in the 
following, and with A — A8 nm, L = 0.34 nm and 
jds =1000 pN. 

• Unless otherwise stated, the double-stranded DNA 
linker is made by N^s — 3120 bps, while the two 
single-stranded linkers are made by Ngs = 40 -f 
n bases each, where n is the number of open DNA 
bases (in other words we included on each side a 40 
bases single-stranded linker). 

• We worked at fixed temperature ksT = 4 pN nm, 
corresponding to T = 16.7 °C. 

• We used the dynamical equations for the polymers 
defined above, Eqs. within the discrete proce- 
dure illustrated in section fill Bl and with transition 
rates (j26|) for the fork with attempt rate r = 10^ 
Hz. 

• The matrices F corresponding to the setups in 
figure |T] are given in Eq. pSj): we used 7^" = 
Iss = 7' = 2 • 10~^ pN s/nm. We used a value 
7 = 1.67 • 10"^ pN s/nm for the viscosity of the 
beads. 

• The time step was fixed to At = 10~* s; this value 
ensures a correct integration of the equation of mo- 
tion in all the regimes discussed below. Even if in 
some cases a larger integration step could be used, 
we decided to keep it fixed in order to be sure that 
discretization biases are not present. 

The values of the spring constants ki and fc2 and of the 
force / in Eq. (|24|) varied in different simulation runs and 
will be specified later. 

The program we used for the numeri- 
cal simulations can be downloaded from 



http : //www . Ipt . ens . f r/"^zamponi[ A user-friendly 
version will be made available as soon as possible. 



D. Details of the numerical simulations 

We performed numerical simulations of the molecular 
constructions depicted in figure lU with the following 
specifications: 



E. Limits of validity of the dynamical model 

Our model of the polymer dynamics suffers from two 
main limitations. 
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First we keep only one collective coordinates for each 
polymer (its extension) associated to the longest relax- 
ation mode. Faster modes are discarded. The approxi- 
mation is justified provided there is no other mode slower 
than the typical sojourn time on a base-pair. From the 
discussion of Section [II C li the number of unzipped base- 
pairs, n, cannot be well above a thousand. 
Another upper limit on n comes from the assumption 
that the force is uniform along the polymer. In principle 
the force is a function of the time t and the location 
y along the polymer, which obeys a diffusion equation 
with microscopic diffusion coefficient ~ {xTs)'^l'^ss 
where x™ is the length of a monomer and t™ = 
its relaxation time. Assume that, at time 0, a base-pair 
closes and the polymer is stretched at the extremity a; = 
by x™. Then the force, initially equal to f{x,t = 0) = 
k™ x™6(x) will decay following the Gaussian diffusion 
kernel. At time t the force density at the extremity is 
f{x,t) — k^Jl x™gl ^I-kD""^], t. The relaxation is over when 
this force excess is of the same order of magnitude as the 
typical thermal fluctuations Sf calculated in ([5]), that is, 
for times 

^ ft., y-^ss) r„ _ 2 10~i° n ps . (35) 
2Tr ksT " f \ J 

When n ~ 1000 the corresponding relaxation time is of 
the order of the sojourn time on a base. 
In conclusion our dynamical model is adapted to ssDNA 
polymers whose length ranges from a few hundreds to 
a few thousands bases. Shorter polymers can be con- 
sidered at equilibrium, while longer polymers cannot be 
modelled without taking account the space-dependence 
of forces. A simple way to tackle this difficulty consists 
in arbitrarily cutting long polymers in 1000-base long seg- 
ments, each modelled as above. This procedure will be 
followed in Section IV Al 

IV. UNZIPPING AT FIXED FORCE 

A. Quasi-equilibrium unzipping 

Before turning to the more interesting case of out-of- 
equilibrium unzipping, we focus on the case of a small 
molecule which is subject to a constant force close to the 
critical force. In this situation the molecule is able to 
visit all the possible configurations. 
We performed a set of numerical simulations at constant 
force / = 16.45 pN, with the setup described in figure[T}3. 
The DNA molecule is a uniform segment of = 500 
base-pairs, with pairing free energy G(n; B) — ngo and 
go = 2.69 fcfiT. The entropic free energy per base of the 
two open single strands is 2gss{f) — 2.684:kBT, there- 
fore the infinite molecule would stay close; we are sligthly 
below the critical force. To the right and left open por- 
tions of the molecule, two single-stranded DNA linkers 
of N^^ = 40 bases each are attached; therefore the total 
length of the single-stranded linkers is Nss = N^g + n, 




FIG. 4: Bottom: Average fraction of time spent on each base. 
The full (blue) curve corresponds to Eq. (|34|) while the dashed 
(black) curve corresponds Eq. (|34|) without the saddle-point 
corrections (the square-root term). The dot-dashed (red) line 
is Peqiji) oc exp[— n Ag] with /S.g — 0.006. 
Top: Effective rates (squares and triangles) estimated from 
the maximization of the probability in Eq. (|36p (r = 10^ 
Hz) without saddle point corrections (full curve of the lower 
panel). The dashed lines are the asymptotic values of the 
rates, see text. We do not report the rates corresponding to 
the full Eq. p4p since they are essentially independent of n. 



where n is as usual the number of open base-pairs. The 
leftmost linker is a double-stranded DNA of Niig = 3120 
base-pairs, whose presence is however irrelevant for the 
scope of this section. The total length of the simulation 
was of T = 7200 s, i.e. of two hours. 



1. A test of the model 

The average fraction of time spent on each base, cor- 
responding to the equilibrium probability distribution 
Peq{n), is reported in the lower panel of figured We 
expect that in the large n limit, Peq{n) ~ exp[— 71(170 — 
2gss{f))] = exp[-nAg], with Ag - 0.006. This is ex- 
pected to break down when Nss is so small that the 
second-order corrections to the saddle-point in Eq. ([M]) 
become important. As it can be seen in figure [H the 
exponential form correctly describe the data. 
We performed additional simulations in which the 



11 



square-root term in Eq. (p4)) was removed. As one 
can see, in this case the small n deviations are much 
more pronounced. It is worth to note that for a non- 
Gaussian polymer one expects a deviation from the ex- 
ponential form at small enough n. However, this analysis 
shows that taking into account the small n corrections to 
W{x, n) systematically reduces this effect. Estimating its 
real order of magnitude therefore requires an exact ex- 
pression oiW{x,n), which could be in principle obtained 
from the recurrence equation l|30p . However this is a 
complicated numerical task that goes beyond the scope 
of this article. What we want to stress here is that the 
inclusion of the square-root term in Eq. (p4)) gives signif- 
icant differences when n < 200 and should therefore be 
included if one wants to analyze the unzipping of small 
molecules. 



B. Out-of-equilibrium opening 

For long molecules the barrier between the closed and 
open states may become very large e.g. ~ 3000 k^T for 
the 50000 bases A-DNA at the critical force /c = 15.5 
pN [31}. The time necessary to cross this barrier is huge, 
and full opening of the molecule never happens during 
experiments. To open a finite fraction of the molecule 
the force has to be chosen to be larger than its critical 
value. The opening can then be modeled as a transient 
random walk, characterized by pauses at local minima of 
the free energy and rapid jumps in between [Tg} . 

1, Analytical calculation of the average time spent by the 
fork on a base 



2. Effective dynamics of the fork 

In a situation where the linkers are short, such that their 
relaxation time is faster than the mean time spent on a 
base, the linkers are able to reach equilibrium before n 
changes. Therefore one might hope to define an effective 
dynamics for the fork, where n changes according to ef- 
fective rates that depend on the variation in free energy 
of the setup on closing or opening a base. 

To this aim we considered the model for the fork dynam- 
ics described in section III C[ but assuming n-dependent 
opening and closing rates. Within this model, the prob- 
ability of a trajectory of the fork is a function of the 
number of upward (M„)/downward (d„) jumps, and the 
time spent on base n, tn'. 



N 

PeffHt)] - l[{rtff{n)Atr-{rlff{n)Atr- 
X (1 - Atirffin) + r^^"{n))Y" 



Given the values of Umdn,tn measured along our tra- 
jectory of duration T, we can infer the effective rates by 
maximizing the above probability. Assuming r'^'^^ At ^ 1 
we obtain 



Gonsider first the case of a fixed force acting on the fork 
while all the other components are at equilibrium as in 
Section III Gl In the transient random walk the opening 
fork spends a finite time around position n before escap- 
ing away and never coming back again in n. The number 
Un of opening transitions n ^ rt + 1, is stochastic and 
varies from experiment to experiment, and base to base. 
The total number of times the fork visits base-pair n be- 
fore escaping is given by the sum of the number u„ of 
transitions from n — 1 to n and of the number Un+i — 1 
of transitions from n + 1 to n. The average time spent in 
n is therefore 



tr. 



(Un) + (Un+l) - 1 

To + rc{n) 



(38) 



where 1/ (ro + rdn)) is the average time spent in n before 
each opening or closing step. Let us introduce the prob- 
ability El^^i of never reaching back position n starting 
from position n + 1. The probability P of the number u„ 
of opening transitions n n+1 during a single unzipping 
simply reads 



PK) = (i-i?;Vi)""^' K+i 



(39) 



From equation (|39p we have that the average number of 
openings of bp n is 



„ ^1 



(40) 



r^ff 



in) 



dn 
in 



(37) 



as estimates for the effective rates. For the full expres- 
sion (|34p. the rates are almost independent of n; on the 
other hand, if the first order correction are neglected, 
one obtains n-dependent rates, consistently with the ob- 
servation that Pf,q{n) is not exponential. These are re- 
ported in the upper panel of figure |3) In both cases, the 
rates are consistent with the detailed balance condition 
rlSf{n)P,q{n) = /(n - l)Peq{n - 1). 



We are thus left with the calculation of E^^^. For infinite 



force £;;^+i 



1 since the fork never moves backward. For 
finite force we write a recursive equation for the probabil- 
ity that the fork never comes back to base n starting 
from base m (> n -I- 1) 



where 



- 9m-E," _i + (1 - qm)El!;^^i 



(41) 



(42) 
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FIG. 5: Top; Average time spent by the fork on position n. 
Bottom: Time spent by the whole setup at extension com- 
prised between xz and xz + Ax, with Aa; = 0.5 nm. 
The black line in both figures represent the theoretical pre- 
dictions from Section riV B II Red points are the results from 
the simulation. Standard deviations are represented by error 
bars in the top panels, and by the thickness of the red curves 
in the bottom panels. 



is the probability to close base n and 1 — (/„ is the proba- 
bility to open it at each step. Note that for forces larger 
than the critical force, we have q-a <\'- the random walk 
is submitted to a forward drift and is transient. The 
boundary conditions for equation ([1T|) are i?" = and 
-E",. = 1 for m oo. 

For an homogeneous sequence the escape probability is 
E — (\ — 2(7)/ (1 — q). For an heterogeneous sequence 
by defining pj^ — ^„'" we obtain the Riccati recursion 
relation 



0; p: 



m+l 



1 — ^m+l 
1 - gm+1 P'^ 



for n > m . (43) 



Equation (j43[) can be solved numerically for a given se- 
quence. The escape probability starting from n -I- 1 is 
then 



- 



n 

m>n+l 



(44) 



and the average time spent in the base n is then obtained 
from (gOl) and 



2. Results from the dynamieal model 

To check whether these theoretical predictions are af- 
fected by dynamical fluctuations of the bead, linkers, and 
unzipped strands we have carried out simulations with 
the model of Section Hill We have carried out 160 unzip- 
pings of the A-phage sequence at the force of 17 pN for 
T = 100 seconds (physical time), with the same molec- 
ular construct of section IIV Al {N^s = 3120 base-pairs 



dsDNA linkers on a side plus 7V£, = 40 bases of ssDNA 
linker at each side of the DNA to be open). For such a 
construct the equilibrium extension of the polymers for 
n open base-pairs is 2Nsslss + N^slds where Ids = 0.3337 
nm,and Iss = 0.4758 nm, and Nss = + n. The stiff- 
ness of the polymers is 1/fce// — Ngs/k"^ + Nds/k'^g with 
A:™ = 160.5 pN/nm and k'Jl = 1450 pN/nm. The relax- 
ation times of the polymers are of the order of 0.1 ms 
for about 400 unzipped bases and 1 ms for about 2500 
open bases, and are larger than the characteristic times 
to open a weak base, of about 2 10~^ s, and a strong 
base, of about 10~^ s. 

We plot in Figure [5] the average time spent by the fork at 
location n for two portions of the sequence, corresponding 
to about 400 and 2500 open base-pairs. The agreement 
between the theoretical and numerical estimates of the 
times is excellent, meaning that the fluctuations of ex- 
tensions of the polymers and the dynamics of the bead 
induce negligible changes on the rates of opening and 
closing, as seen close to the critical force in Section ITVAI 
As experiments do not give direct acces to the time spent 
by the fork at location n we show in Figure [5] (bottom) 
the time t{x3) spent by the unzipped ssDNA between 
extensions xs and X3 + dx. These times are compared to 
their values assuming that the positions X3 of the beads 
are randomly drawn from the equilibrium measure, 



(45) 



where tn is calculated from ([SS]) and P(a;3|n) is calculated 
from an argument similar to those used in section FlIB II 
and can be written up to quadratic order around the 
saddle point as: 



P{x3\n) = 



X e 



2tt 



(46) 



The agreement is, again, excellent. 

Figure O and equation ((15)) show that t{x3) gets con- 
tributions from the times spent by the fork on a set of 
bases whose number depends on the magnitude of the 
equilibrium fluctuations of the linkers. These equilib- 
rium fluctuations increases with the length of ssDNA e.g. 
6x3 ~ 5 nm for 400 unzipped base-pairs and 6x3 ~ 12 
nm for 2500 unzipped bases. Therefore, as the number n 
of unzipped base-pairs increase, the characteristic curve 
of t{x3) gets more and more convoluted (compare left- 
bottom and right-bottom panels in [5]). 
In Figure [6] we compare the value of the ssDNA exten- 
sion from one unzipping, 0:3, to its average value at equi- 
librium, 0:3*, as a function of the number of unzipped 
base-pairs n. The fluctuations in the extension are com- 
patible with the equilibrium deviations. Again no clear 
out-of equilibrium effect is observed. The reason is that, 
even if the single strand is not relaxed in the opening 
time of a base, the fork goes back and forward around a 
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nb. n of unzipped bp nb. n of unzipped bp 



FIG. 6: Total extension of the setup in figure [T)3 at fixed 
number n of unzipped bases for a single unzipping (black line). 
If the fork visits twice or more the same base n we plot the 
average of the extension values. The gray strip represents the 
average value at equibrium, x'i^{n), and the standard devia- 
tion around its value at equilibrium. 

given location before moving away. Therefore the quan- 
tities we have measured are averaged on the number of 
times a base-pair is opened and are close to their mean 
value even in a single unzipping. This can be deduced 
from figure [5] by comparing the total time spent on a 
base (points) with the time to open a base (dashed lines) 

V. UNZIPPING AT FIXED EXTREMITIES 

A. Correlation functions 

One of the main advantages of considering the dynam- 
ics of the linkers and of the beads is that it allows us to 
compute autocorrelation functions and to explore the in- 
teraction between the different parts of the setup, a task 
which would be impossible from a priori calculations. 
We have performed a few simulation with the setup 
shown in figure [T]A_ where the spring constant of the first 
optical trap of extension xi is 0.1 pN/nm and the second 
{xi) has stiffness 0.512 pN/nm. The molecule in the fork 
is uniform with go ~ 2.69fcBT. The only parameter that 
is varied across simulations is the distance between the 
optical traps and thus the typical number of open bases. 
In figure [7] we show two typical cases. What is evident is 
that the single strand has two time scales: one which is 
proper to the fluctuations at n fixed and another which 
is of the same order of magnitude as the correlation time 




FIG. 7: Correlation functions for the setup in figure [TJA. at 
two different values of the number of open bases, Nss = 40-|-n. 



of the fork. As the number of open bases grows the fast 
time scale also grows until it becomes impossible to dis- 
tinguish the two. 

As remarked in section Fill El our model cannot in prin- 
ciple be used when the linkers are made oi n ^ 1000 
monomers. To check for the importance of force propaga- 
tion effects, we ran a simulation for Nss ~ 9700 (bottom 
panel of figure [7|) where we cut each linker in 9 subunits 
of 1000 bases each, plus a final unit which is connected 
to the opening fork. Overall, the correlation functions 
are not much affected by this modification and in partic- 
ular the correlation times are unaffected within numeri- 
cal errors. The main effect of cutting the long linkers is 
that the correlation function of the linker becomes more 
strectched (i.e. if they are fit with exp[— (t/r)'''], the 
exponent (is is sligthly smaller). This is to be expected 
since by cutting the polymer we include more relaxation 
modes, each with its relaxation time. A wider distribu- 
tion of relaxation times implies a smaller exponent (3s. 
In table |V] we compare the results of the numerical simu- 
lation with the predictions of section III B II which do not 
take into account the interactions between different parts 
of the setup. While the simulated result for the single- 
stranded and the double-stranded DNA are not too far off 
from the prediction, the two springs show a much greater 
deviation from the theoretical estimates. This prompted 
us to analyse further the relationship between the fork 
and the bead position as will be discussed later. 
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Theoretical (s) 


Numerical (s) 


Single strand 


4.83 ■ IQ-'^^N^, 


5.4- IQ-^^N^ 


Double strand 


4.96 ■ 10"^ 


~3 - 10"^ 


Spring xi 


1.67- 10"* 


~ 1.5 ■ 10"^ 


Spring X4 


3.26 ■ 10"'' 


~ 7 - 10"^ 


Fork iV,, 


oc 14.2 + 0.013iVss 


1.3- 10"^ +8.4- 10-^iV.,s 



TABLE V: Comparison between the correlation times of the 
setup in figure [TJ\. as computed for an isolated element and 
the result of a complete numerical simulation. In the case of 
the fork, we reported as theoretical value l/kejf, that must 
be multiplied by a viscosity to obtain the relaxation time; it 
turns out that a viscosity ~ 8 • 10~^ pN s/nm matches the 
theoretical and numerical results. 
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le-03 r 

T - 
le-04 r 

le-05 r 




le+04 



FIG. 8: Relaxation times of the correlation functions in figure 
[7] as a function of the number of open bases. In the case of 
the single strand (ss), only the fast relaxation time is plotted. 
For the fork and the single strand, dashed lines indicate a fit 
to T„ = A + BNss (with A = 1.3 ■ 10"=^ and B = 8.4 • 10"^) 
and Tss = CN% (with C = 5.4 • 10"" s). For the others, fuU 
lines are guides to the eye. 
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FIG. 9: Top: Autocorrelation functions for the setup in fig- 
ure [TJ\ when the molecule to unzip is a block-copolymer com- 
posed of alternating stretches of ten strong pairs and ten weak 
pairs. This way the fork correlation time is greatly increased 
allowing us to view effects on the two traps of different optical 
stiffness. 

Botttom: Correlation functions between one of the two beads 
and the number of open base-pairs. Values have been nor- 
malised so that the value at zero time difference is p = 



Mutual information between bead position and 
fork location 



The potential acting on the fork position, in the case of 
an uniform molecule is dictated by the stiffness of the 
rest of the setup only as seen in section III B II That is to 
say that n experiences an harmonic potential with spring 
constant proportional to keff', this in turn predicts cor- 
relation times that are proportional to jp^j which has a 
linear dependence on n. This behavior is in very good 
agreement with the data that has been extracted from 
numerical simulations. 



Figure [Hj shows the dynamical correlations of the fork 
and beads positions. The two beads have different corre- 
lations functions due to the difference in their stiffnesses: 
k ~ 0.5 pN/nm for bead 1 and k ~ 0.1 pN/nm for bead 
2. After an initial decay (taking place over a time propor- 
tional to 1/fc from Section 111 C 3P the bead correlations 
exhibit a quasi-plateau behavior whose height is roughly 
proportional to 1/k. The plateau reflects the correlation 
between the motion of the bead and the one of the fork 
on time-scales of the order of the equilibration time of 
the fork. It appears that soft beads allows one to track 
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the location of the fork better than stiffer beads. 
In the following we will give a closer look to the depen- 
dence of this correlations to the optical trap stiffness, to 
do so we construct a setup as in figure[Tl^, but where the 
stiffness of the optical trap on the left is kept constant 
at 0.512 pN/nm while the stiffness of the one the right is 
varied across two orders of magnitude^. 
To give quantitative support to this statement we define 
the mutual information / between the position of the 
bead in the optical trap, X4, and the number of open 
base-pairs, n: 

/(.4,n)=5:|dx4F(x4,n)log(^^^^il^), (47) 

where P{x4^,n) is the joint probability density for the 
bead to be at position 0:4 while there are n open base- 
pairs; P(n) and P{x4) are the two marginals. Note that 
the definition of mutual information does not suffer from 
the problems which arise with entropy when we switch 
between a continuous and a discrete definition, that is 
to say that binning with sufficiently small bins does not 
change the mutual information. 

/ can be easily computed by keeping track of the times 
passed at a given bead position and given number of open 
bases during a run of the simulation. As stressed before 
the fact that the 0:4 coordinate must be binned has neg- 
ligible effects on the computation of entropy. For very 
large stiffnesses the amplitude of the oscillations of the 
bead can become very small, and thus a lack of sensi- 
tivity in the measure of the position of the bead could 
become an issue. Fortunately current state of the art in 
optical trap cannot attain stiffnesses larger than, say, 1 
pN/nm with micrometer beads [3 33 . In this regime the 
fluctuations of the bead are dominated by the stiffness of 
the trap and thus we can say that {5x\) ~ (/3fc2)~^, see 
Eq. Comparing the fluctuations of the bead posi- 

tion with the sub-nanometer precision A over its location 
yields 

^M^lO- 50, (48) 
which is much larger than unity. 

Figure [10] shows that the mutual information / depends 
only weakly on the sequence but depends strongly on the 
stiffness k of the trap. This behavior can be understood 
very intuitively. Right after a base-pair opens or closes, 
the whole setup in a fixed force experiment has to give 
way; the less an element of the setup is rigid compared 
to the rest, the more it will accomodate for the change 
in n. 



^ The attentive reader might have noticed we changed the stiffness 
of the right bead compared to what it was in the previous section, 
the rationale behind this choice is to keep its value at the center 
of the range in which we will vary the other. 
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k (pN/nm) 

FIG. 10: Mutual information I between X4, and n as a func- 
tion of the trap stiffness, k. Black circles are computed on 
an uniform sequence, while red squares are measured on the 
sawtooth potential described in the caption of figure [9] 



We conclude that, in a single measurement, soft traps 
give more information on the fork location than stiff 
traps. However I is the mutual information between 
the fork and bead locations per measure. As we have 
seen in Section [VAi the correlation times extracted from 
the simulations decrease with k and, as k grows, more 
and more uncorrelated measures can be done in the same 
amount of time. It is thus expected that the information 
per unit of time is not maximal for small values of k. 
In other words, stiffer traps give worse-quality but more 
frequent signals on the location of the fork. Finding the 
optimal value of k would require a detailed analysis of 
the correlation times of the bead and of the fork. In par- 
ticular the size of the bead would affect the optimal value 
for k through the viscosity coefficient, but not the infor- 
mation per measure, /. However this dependence should 
not be crucial since the bead size cannot be varied much 
in experiments: it can be neither too small to exert a suf- 
ficient force, nor too large due to the size of the physical 
setup. 



VI. CONCLUSION 

This paper has been devoted to the presentation of a dy- 
namical model for the different components of the setups 
used in the unzipping of single DNA molecules under a 
mechanical action. Compared to previous studies our 
model does not assume a priori that the polymers in the 
molecular construction are at equilibrium but takes into 
account their relaxation dynamics. It is important to 
stress out that the dynamical description for the linkers 
and the unzipped part of DNA is coarse-grained: the ba- 
sic unity are the polymers themselves, not the monomers 
they are made of. 

As a consequence each polymer is associated with a 
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unique relaxation time. The assumption is justified as 
long as those times are comparable to the typical open- 
ing or closing time of a single base-pair. Longer poly- 
meric chains e.g. ssDNA strands with a few thousands 
bases need to be modeled in a more detailed way; more 
precisely, they should be divided into short enough seg- 
ments along which the force can be considered as uniform 
on the time scales associated to the fork motion. Al- 
though in this paper we did not observe any important 
force propagation effect, these might be more important 
in strongly nonequilibrium situations such as opening at 
constant (high) velocity. We plan to simulate unzippings 
with such molecular constructions in the next future to 
understand how force propagation across the polymeric 
segments can affect the effective rates for closing base- 
pairs in such situations. 

One of our results is that one has to be very careful 
with the expression of the free-energies (entering the dy- 
namical rates) for short polymers, be they linkers or 
ssDNA unzipped strands. Use of the free-energy per 
monomer, obtained from force-extension measures on 
long molecules, as usually done in the literature can lead 
to erroneous results. We have shown that finite-size cor- 
rections to the energetic contributions and the dynamical 
rates have to be taken into account. 
As a main advantage the code we have developed is ver- 
satile: we can easily change setups, for example use fixed- 
force or fixed-position ensemble, and change the number 
and types of linkers and of traps for the beads. We have 
found that, in fixed-force unzippings, the opening and 
closing rates for the fork are not affected by the force 
fluctuations coming from the polymeric chains. For small 
linkers and number of unzipped base-pairs, indeed, force 
fluctuations are large but fast, and are averaged out on 
the characteristic opening-closing time of a base-pair. For 
large linkers or number of unzipped bases force fluctua- 
tions are slow but small, and therefore do not change the 
dynamic of the opening fork. We have also performed 
unzipping simulations at large forces where the opening 
dynamics is transient, and found that the average time 
spent by the unzipped strands at a given extension is 
accurately predicted from the time spent by the fork on 
a base convoluted by the equilibrium fluctuations of ss- 
DNA. Moreover the extension between the extremities 
at a fixed number of open base-pairs in a single unzip- 
ping experiments is compatible with equilibrium fluctua- 
tions of ssDNA and linkers. The program could be easily 
adapted to unzipping at constant velocity, where non- 
equilibrium effects are likely to be more important. 
Our study suggests that one measure of the position of 
the bead in soft traps gives more information on the 
location of the fork than in the case of stiffer traps. 
This statement is however to be considered with caution. 
Beads in stiffer traps reach equilibrium on shorter time 
scales, and the overall rate of information per unit time 
could be higher in stiffer traps. While purely qualitative 
at this stage, such a statement is relevant to the study 
of the inverse problem of unzipping, that is, inferring 



the sequence of the DNA molecule from the unzipping 
signal. We hope that the present dynamical modeling 
will be useful to assess the rate at which information on 
the sequence could be acquired from mechanical single 
molecule experiments. 
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APPENDIX A: LANGEVIN DYNAMICS OF 
COUPLED POLYMERS 

One of the simplest models of polymer dynamics is the 
one proposed by Rouse [44| . where the polymer is de- 
scribed as a chain of beads which are modeled as Brow- 
nian particles, linked by harmonic springs. 
While it is true that this model is very crude because 
it ignores hydrodynamic interactions and excluded vol- 
ume effects, it has the huge advantage of being largely 
solvable. Therefore we will use it now as the basis for a 
few considerations that will be then generalized to more 
realistic models. 

Our aim is to write a system of coupled equations for the 
time evolution of a certain number of marked points on a 
(hetcro)polymcr. One of these points will be for instance 
the location of the opening fork. In the case of a double 
DNA strand attached to a single strand, one point will 
mark the location where the two different polymers are 
attached (see the examples in figure [T]). Note that if the 
marked points we focus on are far apart, only the slower 
modes of the system will be relevant, as the fast modes 
describe local relaxations of the chain. Therefore in the 
following we want to focus on a long wavelenght/long 
time effective description of the chain. 



1. The dynamics of a single polymer 

a. The model and its normal modes 

As the simplest case we consider a polymer composed 
of N identical springs, each with an identical link at one 
end. The first is connected to a wall that has infinite mass 
(or, better still in this framework, infinite viscosity) and 
on the last is exerted a force /. The Langevin equations 
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describing such a polymer can be written as: 

7„iUi = -2fc,„Mi + kmU2 + m 
* 7mUn = ^SfcjnUn + fcrn^^n-l + ^mW„+l + Tj^ (-A-1) 



.ImUN = -k„iUN + kmUN-1 + f + VN , 

where the r]i are white Gaussian noises of zero mean and 
variance 



(A2) 



Let us for the moment neglect the noise term. Then, 
defining = Jm/km, we can rewrite formally these 
equations as 

TmUn = -2u„ + U„-l + U„+l , Vrt , (A3) 

supplemented by the boundary conditions 

Wo = , UN+1= UN + f /hn . (A4) 

A standard way to find the normal modes of the linear 
system above is to search for solutions of the form Un (t) = 
M„(0) exp(— Zt/rm). One can easily show that the general 
solution satisfying the first boundary condition mq = 
has the form 



Un{t) cx svLv{qn)exp{-l{q)t/Tjn) , 
/(q) = 2(l-cos(g)) . 



(A5) 



The second boundary condition (jA4p requires that 
UN+iit) — — f /km — const. Since we can always 

add the constant value to UN+i{t)^ we can replace this 
boundary condition by UM+i{t) — UN{t). This requires 
that sin(g7V) ~ sin(q(iV + 1)), then q = [-k /2 + pi:) / N . 
The slowest mode then correspond to q = 7r/2/iV, that 
for large TV gives a relaxation time 

t{N) = T„Jl{^/2/N) ^ -^TmN^ . (A6) 

which proves the validity of the scaling in Eq.ipO]). 



b. Recurrence equations for a fixed end 

We want now to write a system of coupled equations for 
a certain number of points on the polymer by integrating 
out the u's we are not interested in. To begin, we focus 
on the end point un. 

It is convenient to perform a Laplace transformation and 
write 



/•OO 

Jo 



Then Eq. (|A5p becomes in Laplace space 

(2 - l)Unil) = U„+l(/) + U„-l(0 , 



(A7) 



(A8) 



with the same boundary conditions uq(1) = 0, and 
un+i{1) — ujv(0 — [f /Kn)5{l)- For / ^ the latter 
condition reduces to un+i{1) = un{1) as discussed above 
for the normal mode analysis. 
We introduce a function 

C«-l(0 ="n-l(0/"«(0 • (A9) 

Substituting the latter relation in (|A8[) we get 

(2 - / - Cn-l(0)"n(0 = W„+l(0 , (AlO) 

from which we get a Riccati recurrence equation 



Co(0=0 (due to Mo = 0) 



(All) 



This recurrence can be solved and the function Cn(0 com- 
puted for all n. 

Since we are interested in the large time limit, we can 
expand the function Cn(0 fc)r small I; we obtain 



Cn(0 = 



n n(l + 2n) 



n + l 6(1 -l-n) 
n(6 + 19n+16»^ + 4n3) ,, 
180(1 + n) 



(A12) 



One obtains the effective equation for by substitut- 
ing the expression above in (jAlOp and setting n = N. 
Keeping only the linear term in I and the leading terms 
in ^ 1, we get 



1 N\ 
1 + - ) un{1) = UN+i{l) 



Moving back to time domain, we obtain 



N . 



1 

— ujv + (ujv+i - um) , 



(A13) 



(A14) 



which is equivalent, using the boundary condition ujv+i — 

UN = f/km, to 



(A15) 



In this way we got an effective equation for the endpoint 
of the polymer that is still a linear first order differential 
equation and takes into account only the slowest mode 
of the chain. 

There is however an inconvenient: in fact a straightfor- 
ward computation shows that the relaxation time ob- 
tained from Eq. (|XT5)) is t{N) = t^N^/S, that differs by 
a factor 7r^/12 from the correct value given by Eq. (jA6[) . 
The origin of this discrepancy is clearly in the fact that 
the expansion we made in Eq. (|A12[) is not convergent at 
fixed / for n ^ cxd, as successive terms in the series are 
of order ti^P'^Ip. 

Let us then go back to the computation of the normal 
modes of the system within this formalism. The second 



18 



Xi 

I I I 



Ui 



I I I I- --M- 

Uff-l l'-(J\/-l) V-2 V2 Vm-1 



FIG. 11: Two joint polymers subjected to an ex- 
ternal force /. xi marks the endpoint of the 
first polymer made of A'' links whose endpoints are 
ui,U2,-- - ,un-i,un = x\. The second polymer originates 
from xi and is made by 2M — 1 links, whose endpoints are 

U_(A./_l),W_(M-2),- ■ ■ ,V-l,Vl,-- - ,VM-2,VM-1,X2. 



boundary condition UN+iil) = WAr(0 implies Cn{1) = 1- 
The normal modes are the solutions of this equation with 
respect to I. One can show from the exact expression of 
CnIi) that 



lim NiCNiq'/N^) - I] = -q cot{q) ee ({q) . (A16) 



The zeroes of this function are q — 7r/2 + kn; there- 
fore the solutions of CwO = 1 tend for large N to 
I — {tt/2 + pir)'^ /N'^, in agreement with the exact re- 
sult of the previous section. An inspection of Eqs. (jA12p 
and (|A16[) shows that the small I expansion of Cn{1) is 
equivalent to perform a small g expansion of C.{q) in order 
to find its first zero. This indeed yields ((q) ~ —1 +'q^/3 
that gives q — \/3 for the first zero that gives back 

t{n) = TmNys. 

Then one can check that a higher order expansion in I 
(or equivalently in q) produces a more accurate result; 

indeed the series of C(5) converges for g < tt while the 
zero is located in 9 = 7r/2. It is easy to show that if one 
truncates the series to order p, the difference between the 
solution and the true zero is exponentially small in p. 



Discussion 



The conclusion of this section is that Eq. IIA15\) is a cor- 
rect description of the dynamics of the end of the polymer 
in the limit of large N and large times. While it captures 
the correct scaling with N of the relaxation time, the co- 
efficient is wrong by a factor of 7r^/12 ^ 0.82. Still this is 
quite satisfactory for our purposes since the experimental 
error in the determination of is of the same order of 
magnitude. Better approximations can be obtained by 
truncating the expansion of Cn{1) to higher orders in 
therefore obtaining a higher order differential equation 
for UN(t). 

In the following, we will derive the coupled equation for 
many marked points along the chain, limiting ourselves 
to the first order truncation. This produces first-order 
differential equations of the Langevin type. 



2. Dynamics of two coupled polymers 

We will now show how to use this formalism to derive cou- 
pled equations for different points on a composite poly- 
mer. We keep neglecting the noise, that we will reintro- 
duce at the end of this section. 

As a simple example, let us consider the polymer drawn 
in figure [TTl It is composed by N monomers of type "U" 
linked to 2M — 1 monomers of type "V" . The two types 
of monomers might differ for the value of the microscopic 
sping constant, bead viscosity etc. If the monomers are 
identicals, then we are just marking a point in the middle 
of a polymer. 

The effective equation for the endpoint of polymer U can 
be derived following the analysis of the previous section. 
We denote xi = upf and we get 

N 

(A17) 

where the last term is the "external" force that the poly- 
mer V exerts on U. 



a. Integration of the V polymer 

Now we want to integrate out all the monomers 
• • • ,um-i in order to obtain the coupling be- 
tween xi and X2- To this aim, and in order to keep 
the formalism symmetric, we can start from the middle 
of the polymer V by integrating simultaneously w_i and 
vi in order to obtain effective equations for v-2 and V2, 
and so on. In Laplace space (note that now in Eq. (|A7p 
Tjn = T^) the equations for v±i have the form 



(2-0«-i(0 = «-2(0 + ^^i(0 , 

i2~l)vi{l)^V2{l)+V-l{l) . 



(A18) 



These can be easily solved to get v±i as function of v±2- 
Iteration leads to the following form for the equation after 
n steps: 



Cn(0 V-n-lil) = W-«-2(0 + Vn{l) Vn+l{l) 
^n(0 Vn+l{l) = Vn+2{1) + Vn{l) V-n-l{l) ■ 



(A19) 



One can check that this form is stable under one step 
of iteration and the following recursion relations are ob- 
tained: 



^0 = 2 
% = 1 

e 



ri+l 



2-1 



(A20) 



Vn+l = JY^. 



where the initial values are determined by consistency 
between (jAl8P and (|A19P for n = 0. These recurrences 
are easily solved by introducing the two quantities An = 
^/(^.n — Vn) and B„ = 1 / (^„ -|- ?]„ ) ; these satisfy the same 
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recurrence in (jAlip except for the initial condition which 

is different and determined according to (|A20|) . 

At the leading order in n ^ oo and at first order in /, we 

get 



Vn{l) = 



1 



2n, 



- —I 



-I 



(A21) 



1 

2n 
1 2n, 

Finally one obtains from this procedure a coupled equa- 
tion for and vm-i where also xi = V-m and 
X2 = vm appear. 

b. Coupled effective equations 

To obtain the coupled effective equations one starts from 
the following system: 



Im 3 



Xi 

£.M-2{1) V-M+l{l) = Xi 

^M-2(0 VM-l{l) ^X2+ riM-2{l) V-M+l{l) , 
(1 - 1)X2 = VM-1 + f , 

(A22) 

where the first equation is just the Laplace transform 
of Eq. (|A17|1 (recall that we are using the definition of 
Laplace transform (|A7|1 with r,„ = t^), the second and 
third equations are Eq. (jA19[) for n = M — 2 and the last 
equation is the Laplace transform of the equation for X2, 
that in time domain reads ^^1X2 = —k}^{x2 — vm-i) + f ■ 
Eliminating v-m+i and vm-i from these equations, us- 
ing the recurrence equations (|A20p and the result (jA21[) 
we finally get the coupled equations: 



N 



■x\ + 



riM-2{l) V-M- 



Xl) 



(^UN,V2M^ 



V 2M ■ I V 2M ■ 



•^'1 6 "^2 - 

2M V^2 



(A23) 

At this point we reintroduce the free energy of the poly- 
mer chain, defining Ni = N and N2 = 2M - 1 2M : 



F{X1,X2) = ^^x 



(A24) 



3. Beads 

At this point we should add the beads that are used for 
the optical manipulation of polymers. These beads are 
optically tweezed or subjected to magnetic fields in order 
to apply forces on the polymers. In the former case, 
the force acting on the bead is a harmonic force / — 
—k{x — X), while in the latter it is constant, / = fext- 
Each bead is characterized by a friction coefficient that 
can be computed using the Stokes law; we denote it by 
7. Typically they are of the order of 10~^ pN s/nm, 
i.e. much bigger than the microscopic viscosity of the 
polymers 7^ ~ 10~® pN s/nm. 

In presence of a bead attached to the endpoint of a poly- 
mer, the equations of motion (|Aip . (jA18|) . etc. remain 
valid, but one should add to the coordinate describing 
the position of the bead the contribution 7 to the vis- 
cosity. For instance, if there is a bead attached to the 
endpoint un, the last Eq. (jAl|) reads 

(7 + lni)uN = -kmUN + fcrn^W-l + f + VN ■ (A28) 

Then the above derivation still holds because the last 
equation is not used until the end. The only modification 
will be the inclusion of 7 on the diagonal element 
corresponding to the coordinate of the bead. 
Therefore to describe the beads attached to the end of 
the molecular construction in figure [Tl we modify the 
matrix F as above; and in case A, we add to the free 
energy a term ^k(x4 — A")^, while in case B we add a 
term -fextX3- 

In the case of figure [T]A_, one also has to include the left 
bead. In this case, if we call V the first polymer after the 
bead, we can start from a system of equations identical 
to (jA22[) . but with the first equation replaced by 



72:1 = —kxi + k'!^{v- 



M+l 



Xl) 



This will lead again to (|A26|) with 
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(A29) 



(A30) 



and a matrix 



and 



F = 



— I hn 3 ' 'm 3 'm Q 
^VN2 ^VN2 
Im e Im 3 , 



SO that we can write the system above as 



0X1 



(A25) 



(A26) 



where / = (0, /) is the external force vector and we rein- 
troduced the noise term ff that we neglected before. 
The correlation function of the noise at this point is 
determined by the requirement that the fluctuation- 
dissipation relation is verified. This imposes that 



(A27) 



Fix,,X2) = -xl + ^{x2-x,)\ (A31) 



4. Description of a generic setup 

The arguments of the previous section suggest that in 
the general case a bead can be treated "as a partic- 
ular instance of a polymer". In other words, we can 
consider the setups in figure [1] as chains of p joint ele- 
ments U = Ui,U2, ■ ■ ■ ,Up; each element can be an "op- 
tical trap" (i.e. a spring) or a polymer of Ni,N2, • • • , Np 
monomers respectively (in the case of an optical trap we 
set by default Ni = 1). The endpoint of each element is 
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denoted by Xi and x = {xi,X2, • • • , Xp) is the state vector 
of the system (we also define xq = 0). 
The total free energy is then F{x) = J2^=i {xi—Xi^i) 
where VFjy. (.t) = ^kx^ for an optical trap of stiffness k. 
Then Eq. ((X26|) holds, with i , j running from f to p and 
the noise correlation matrix is given by ()A27|) . 
The matrix F must be constructed as follows. Each di- 
agonal term Ta^ related to Xi, is the sum of a Stokes 
term coming from a bead possibly attached to Xi, and 
the contribution coming from the two elements adjacent 
to Xi (except for i = p when there is only one contribu- 
tion): 



-7r 



(A32) 



(the first term is present only if there is a bead attached 
to Xi). All the off-diagonal elements are zero except the 
ones adjacent to the diagonal (i.e. connecting Xi and 
Xi±i) which get a contribution from the polymer con- 
necting these two ends: 



= r 



7i 



6 



1 



(A33) 

Note that this final formulation is independent of the 
Gaussian form of F{x) that we assumed in the derivation, 
therefore we will use it also for non-Gaussian polymers 
substituting the appropriate form of Fix) in Eq. (|A26p . 
To conclude this section, note that a further check of the 
quality of the first order approximation can be done as 
follows. If we consider a single polymer made of A'^i -I- 
bases, the corresponding relaxation time is predicted to 
be T = Tm{Ni + N2)'^/3. On the other hand, we could 
consider two coupled polymers of A^i and N2 bases fol- 
lowing Eq. (|J23|) for yl^^ = 7™ and fc^'^ = fc„. The 
coupled equation can be exactly solved and yield two dis- 
tinct relaxation times (that typically differ by a factor of 
10); the slowest relaxation time can be compared with 
T = Tm{Ni -I- A^2)^/3. We found that the difference is at 
most 20%, and the error is maximal for Ni ~ A''2 while it 
decreases when one of the two polymers is much longer 
than the other. 



APPENDIX B: TRANSITION RATES FOR THE 
FORK DYNAMICS 

We now consider a fork n attached to the polymers. For 
simplicty we consider the case of a single polymer whose 
extension is x and free energy is W{x,n). We want to 
construct a stochastic process that samples the equi- 



librium distribution Peq{x,n) 



-PW{x,n)-G{n;B) 



/z, 



where —G{n; B) is the free energy gain in closing the 
first n bases of DNA, as defined in Eq. ([5]). 
The random process is constructed as follows. The 
Langevin equation discussed in the previous section is 
discretized with time step At. If at a given time t the 
system is in a state (a;,n), we allow three possible tran- 
sitions: 



• {x, n) (a; + Ax, n) with rate H^(x, n, Ax), 

• {x, n) (x + Ax, n + 1) with rate H°{x, n. Ax), 

• (x, n) ^ [x + Ax, n — 1) with rate H^{x, n, Ax). 
We must have 

J dAxH^ix, n. Ax) + H°{x, n, Ax) + H%x, n. Ax) = 1 

(Bl) 

Moreover we can define rates r'*'°'^(a;, n) — 
j dAxH^'°''^{x,n,Ax) that represent the rates to 
stay, open or close n independently of Ax. In a practical 
implementation we first decide whether to open, close or 
stay according to and then extract Ax from the 

distribution H'''°'''{x, n, Ax)/r'*'°'^(x, n). 
The detailed balance conditions read 

P(n,x)H°{x,n,Ax) = 

P{n + l,x -f Ax)H''{n + l,x + Ax, -Ax) 
P{n,x)H''{x,n,Ax) = 

Pin ~ l,x + Ax)H°{n - 1, x + Ax, -Ax) 
P(n,x)iJ'*(x,7i, Ax) = 

Pin, X + Ax)!!" in, x -I- Ax, - Ax) 

(B2) 

We assume that the rate for opening is given by the prod- 
uct of a term that only depends on the binding free energy 
as in Eq. (j2ip and a term corresponding to a standard 
Langevin step: 



H°(x,n,Ax) = ^Ate^("^^)-«("+i^^) 

In 



UnTAt 
X 1 / exp 

V 7« 



ATAt 



Ax - 



f{x,n)At 



In 



(B3) 



Note that integrating over Ax we find r°{x,n) = 
^A^gG(«;B)-G(n+i;S) ^ g-^o(f)„+l.f),^+2) ^ cousisteutly 

with Eq. (PT|) . 

Now it is easy to show that the following expression for 
H'^{x,n, Ax) follows from the second detailed balance 
condition: 

H'^ix, n. Ax) = rAt e0^(-:n)-l3W{x+Ax,n~i) 



UirTAt 
X 1 / exp 

V 7n-l 



ATAt 



Aa 



f{x + Ax,n-l)Aty 

7"-i / 
(B4) 



and that the first condition is then automatically satis- 
fied. Up to now we did not specify the form for f{x,n). 
However for a generic /(x, n) the rate above is not Gaus- 
sian. To obtain a Gaussian rate we assume that 



f{x,n) 



dW{x, n) 
dx 



(B5) 
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and perform the following simplifications assuming that 
At is small: 

y l3W{x,n-l)-f3W{x+Ax,n-l)~Pfin-l,x+Ax) 



Finally, the rate at constant n is simply given by 

H^{x, n, Ax) = [1 — i'°(x, n) — n)] 

f{x,n)Af' 



7«-i 
'4TAi 
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X A / exp 
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7n-l 



U-KTAt 
X ^ / exp 
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+ 2 



/3a2vi7(x,n- 1) 



Ax' 



f{x + Ax,n- l)At 

7n-l 



(B6) 



Neglecting 0{Ax^) one obtains a Gaussian distribution 
for Ax, and computing first and second moment of the 
Gaussian one can see that at the lowest order in At it is 
equivalent to 



H'^ix, n, Ax) = rAt ef^w{x.n)-i3W(x,n-i) 

f{x,n 



UnTAt 

X i / exp 

V 7«-i 



7»-i 
'4TAt 



Aa 



l)At 



7«-i 



(B7) 



From the above expression we deduce that the rate for 
closing is r'=(x,n) = ^At e'^^^^'")"''^'^-"-!); and one has 
first to close, and then perform a Langevin step with force 
f{x,n— 1) and friction 7„-i. 



and it is easy to see that this verifies the third detailed 
balance equation if Eq. (jB5P holds and higher orders in 
At are neglected. 

To resume, the implementation of the algorithm is the 
following: 

1. Choose if stay, open or close, with rates r^'°''^{x, n) 
respectively. 

2. If open, first perform a Langevin step at n and then 
increase n by one. 

3. If close, first decrease n by one and then perform a 
Langevin step at n — 1. 

4. If stay, just perform a Langevin step at n. 

5. Goto 1. 

The extension of the above derivation to a case where 
many polymers are present is straightforward, since the 
only polymers whose rate are coupled with n are the two 
adjacent ones. All the other polymers are not influenced 
by n and one can use standard discretized Langevin dy- 
namics. 
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